## cargar paquetes
import geopandas as gpd
import rasterio
import xarray as xr
import rioxarray
import leafmap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as npAnálisis geospacial con Python
Introducción
Paquetes de análisis espacial:
Geopandas: paquete para datos vectoriales (equivalente a
sf)Xarray: paquete para datos multidimensionales (equivalente a
stars?)Rasterio: interfaz a GDAL para datos raster (equivalente a
terra, pero de bajo nivel)Rioxarray: puente entre Xarray y Rasterio.
Paquetes de visualización interactiva:
Leafmap: interfaz sencilla y unificada para mapas web de múltiples backends (folium, ipyleaflet, maplibre).
MapLibre: mapas 2D y 3D avanzados.
HyperCoast: visualización 3D de datos hiperespectrales con pocas líneas de código.
Paquetes de análisis especializado:
WhiteboxTools: herramientas de geoprocesamiento de análisis hidrológico, análisis del terreno, y análisis de datos LiDAR.
Geemap: puente a Google Earth Engine.
DuckDB: análisis de datos espaciales en base de datos de alto rendimiento.
Apache Sedona: computación distribuida para procesado de datos geoespaciales.
Desarrollo de aplicaciones:
Voila y Solara: creación de aplicaciones con Python y Jupyter.
Shiny: creación de aplicaciones y dashboards.
Vamos a comprobar que los paquetes se cargan correctamente y a crear un primer mapa:
## crear mapa
map = leafmap.Map(center = [40, -100], zoom = 5, height = '500px')
## añadir mapas base
map.add_basemap('OpenTopoMap')
map.add_basemap('USGS.Imagery')
## mostrar mapa
mapAnálisis vectorial con GeoPandas
Vamos aprender:
GeoDataFrameyGeoSeriesLeer/Exportar datos vectoriales
Operaciones comunes
Visualizar datos
Trabajar con CRS
Geopandas
Una diferencia con sf, es que un GeoDataFrame puede tener varias geometrías, pero solamente una está activa en cada momento. La geometría se accede a través del atributo .geometry. Vamos a crear un GeoDataFrame:
## datos
data = {
'City': ['Tokyo', 'New York', 'London', 'Paris'],
"Latitude": [35.6895, 40.7128, 51.5074, 48.8566],
"Longitude": [139.6917, -74.0060, -0.1278, 2.3522]
}
## convertir a DataFrame
df = pd.DataFrame(data)
## convertir a GeoDataFrame
gdf = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df.Longitude, df.Latitude))
## imprimir
print(gdf) City Latitude Longitude geometry
0 Tokyo 35.6895 139.6917 POINT (139.6917 35.6895)
1 New York 40.7128 -74.0060 POINT (-74.006 40.7128)
2 London 51.5074 -0.1278 POINT (-0.1278 51.5074)
3 Paris 48.8566 2.3522 POINT (2.3522 48.8566)
Leer datos:
## url a geojson
url = "https://github.com/opengeos/datasets/releases/download/vector/nybb.geojson"
## leer
data_geojson = gpd.read_file(url)Exportar datos:
gdf.to_file('00_data\\test-data.geojson')
gdf.to_file('00_data\\test-data.gpkg', layer = "layer1")
gpd.list_layers('00_data\\test-data.gpkg')| name | geometry_type | |
|---|---|---|
| 0 | layer1 | Point |
Sistemas de Referencia de Coordenadas
Comprobar CRS:
print(data_geojson.crs)EPSG:2263
Reproyectar o transformar CRS:
gdf_4326 = data_geojson.to_crs('EPSG:4326')Medidas
Para calular medidas como áreas y distancias, necesitamos un CRS que esté en unidades como metros (p. ej. EPSG:3857).
## reproyectar
gdf_3857 = gdf_4326.to_crs('EPSG:3857')
## crear índice
gdf_3857.set_index('BoroName', inplace = True)Calcular área:
## área
gdf_3857["area"] = gdf_3857.area
## perímetro en km
gdf_3857["perim"] = gdf_3857.length / 1000Extraer características geométricas
## centroide
gdf_3857['centroid'] = gdf_3857.centroid
## bordes
gdf_3857['boundary'] = gdf_3857.boundaryDistancias: vamos a ver la distancia desde el centroide de Manhattan al resto de centroides
## centroide de referencia
manhattan_centroid = gdf_3857.loc['Manhattan', 'centroid']
## calcular distancia
gdf_3857['dist_to_manhattan'] = gdf_3857['centroid'].distance(manhattan_centroid)Crear buffers:
## crear buffer (metros)
gdf_3857['buffer'] = gdf_3857.buffer(3000)
## visualizar
fig, ax = plt.subplots(figsize=(10, 6))
gdf_3857["buffer"].plot(
ax = ax,
alpha = 0.3,
color = "orange",
edgecolor = "red",
linewidth = 1,
label = "3km Buffer Zone",
)
gdf_3857.plot(
ax = ax,
color = 'lightblue',
linewidth = 1,
edgecolor = 'navy',
label = 'Original Boundaries'
)
plt.title('Distritos y Buffer')
plt.legend()
plt.axis('off')(np.float64(-8272486.092246463),
np.float64(-8197855.319316324),
np.float64(4931921.070848358),
np.float64(5006269.877978747))
Crear convex hulls:
gdf_3857['chull'] = gdf_3857.convex_hullVisualización
Vamos a especificar parámetros para mejorar las visualizaciones:
plt.rcParams['figure.dpi'] = 150Mapa temático:
## crear lienzo
fig, ax = plt.subplots(figsize = (10, 6))
## añadir mapa
gdf_3857.plot(
column = 'area',
ax = ax,
legend = True,
cmap = 'YlOrRd',
edgecolor = 'black',
linewidth = .5
)
## añadir elementos al mapa
plt.title('Distritos de NYC por área', fontsize = 16, fontweight = 'bold')
plt.axis('off')
plt.tight_layout()
plt.show()Mapa interactivo con Folium:
gdf_3857.explore(
column = 'area',
cmap = 'YlOrRd',
tooltip = ['area', 'dist_to_manhattan'],
popup = True,
legend = True
)Relaciones espaciales
Las relaciones espaciales en Python se estrucutran: objeto1.relación(objeto2). Los métodos intersects, touches, crosses, contains, contains_properly… devuelven un Pandas Series de True/False.
Intersección con Manhattan:
## extraer Manhattan
manhattan_geom = gdf_3857.loc['Manhattan', 'geometry']
## intersección
intersects_manhattan = gdf_3857.intersects(manhattan_geom)
## filtrar los que intersecan
gdf_3857.loc[intersects_manhattan, :]| BoroCode | Shape_Leng | Shape_Area | geometry | area | perim | centroid | boundary | dist_to_manhattan | buffer | chull | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| BoroName | |||||||||||
| Queens | 4 | 896344.047763 | 3.045213e+09 | MULTIPOLYGON (((-8219461.955 4952778.645, -821... | 4.928308e+08 | 360.377573 | POINT (-8217436.761 4969318.623) | MULTILINESTRING ((-8219461.955 4952778.645, -8... | 19456.427458 | POLYGON ((-8229469.036 4965425.107, -8229556.9... | POLYGON ((-8231045.21 4944993.298, -8233479.86... |
| Brooklyn | 3 | 741080.523166 | 1.937479e+09 | MULTIPOLYGON (((-8222843.701 4950893.717, -822... | 3.129684e+08 | 297.784594 | POINT (-8231817.476 4960085.209) | MULTILINESTRING ((-8222843.701 4950893.717, -8... | 19586.763567 | POLYGON ((-8245284.602 4956687.136, -8245296.1... | POLYGON ((-8235790.052 4949053.259, -8237874.4... |
| Manhattan | 1 | 359299.096471 | 6.364715e+08 | MULTIPOLYGON (((-8238858.852 4965914.967, -823... | 1.032201e+08 | 144.706085 | POINT (-8233984.776 4979551.696) | MULTILINESTRING ((-8238858.852 4965914.967, -8... | 0.000000 | POLYGON ((-8237068.209 4963506.169, -8237106.8... | POLYGON ((-8240208.861 4965683.834, -8242627.6... |
| Bronx | 2 | 464392.991824 | 1.186925e+09 | MULTIPOLYGON (((-8226155.114 4982269.852, -822... | 1.929250e+08 | 187.221157 | POINT (-8222783.625 4990631.161) | MULTILINESTRING ((-8226155.114 4982269.852, -8... | 15755.010025 | POLYGON ((-8233207.11 4988220.645, -8233189.61... | POLYGON ((-8224095.48 4980733.074, -8225054.51... |
Tocan a Manhattan:
gdf_3857.loc[gdf_3857.touches(manhattan_geom), :]| BoroCode | Shape_Leng | Shape_Area | geometry | area | perim | centroid | boundary | dist_to_manhattan | buffer | chull | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| BoroName | |||||||||||
| Queens | 4 | 896344.047763 | 3.045213e+09 | MULTIPOLYGON (((-8219461.955 4952778.645, -821... | 4.928308e+08 | 360.377573 | POINT (-8217436.761 4969318.623) | MULTILINESTRING ((-8219461.955 4952778.645, -8... | 19456.427458 | POLYGON ((-8229469.036 4965425.107, -8229556.9... | POLYGON ((-8231045.21 4944993.298, -8233479.86... |
| Brooklyn | 3 | 741080.523166 | 1.937479e+09 | MULTIPOLYGON (((-8222843.701 4950893.717, -822... | 3.129684e+08 | 297.784594 | POINT (-8231817.476 4960085.209) | MULTILINESTRING ((-8222843.701 4950893.717, -8... | 19586.763567 | POLYGON ((-8245284.602 4956687.136, -8245296.1... | POLYGON ((-8235790.052 4949053.259, -8237874.4... |
| Bronx | 2 | 464392.991824 | 1.186925e+09 | MULTIPOLYGON (((-8226155.114 4982269.852, -822... | 1.929250e+08 | 187.221157 | POINT (-8222783.625 4990631.161) | MULTILINESTRING ((-8226155.114 4982269.852, -8... | 15755.010025 | POLYGON ((-8233207.11 4988220.645, -8233189.61... | POLYGON ((-8224095.48 4980733.074, -8225054.51... |
Comprobar que los centroides están dentro de su polígono:
gdf_3857['centroid'].within(gdf_3857.geometry)BoroName
Staten Island True
Queens True
Brooklyn True
Manhattan True
Bronx True
dtype: bool
Relaciones binarias (dos geometrías)
gdf_3857['intersection'] = gdf_3857.intersection(manhattan_geom)
gdf_3857['difference'] = gdf_3857.difference(manhattan_geom)
gdf_3857['intersection'].explore()Verificaciones
gdf.is_empty
gdf.is_valid
gdf.is_simple0 True
1 True
2 True
3 True
dtype: bool